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Abstract. - A density functional theory of two-dimensional freezing is presented for a soft 
interaction potential that scales as inverse cube of particle distance. This repulsive potential 
between parallel, induced dipoles is realized for paramagnetic colloids on an interface, which 
are additionally exposed to an external magnetic field. An extended modified weighted den- 
sity approximation which includes correct triplet correlations in the liquid state is used. The 
£_) ' theoretical prediction of the freezing transition is in good agreement with experimental and 

simulation data. 



> 

A microscopic theory of freezing and melting is a great challenge in statistical physics. 
There are two complementary approaches to the liquid-to-solid transition: first, classical den- 
sity functional theory [1-3] starts from liquid state and views the solid as a condensation of 
liquid density modes, hence it is a liquid-based approach. Second, crystal elasticity theory [4] 
is a solid-based theory where the liquid is viewed as a solid with an accumulation of defects. 
In three dimensions, the freezing transition is first order and it is known that it is not de- 
fect mediated. Here, density functional theory provides a molecular theory for the freezing 
transition. Crystal elasticity theory is appropriate to two dimensions and predicts a possible 
— JL , scenario of two-stage melting via an intermediate hexatic phase [5-8] . The advantage of den- 
sity functional theory is that it can be used to calculate the structure of the solid, whereas it 
is not possible to extract the structure of the fluid out of crystal elasticity theory. 
O ■ An excellent realization of a two-dimensional system is provided by paramagnetic colloidal 

^ , particles in a pendant water droplet, which are confined to the air- water interface [9]. If 
an external magnetic field is applied perpendicular to the interface, a magnetic moment is 
induced in the particles resulting into a tunable, mutual dipolar repulsion between them. 
The corresponding interaction pair potential u(r) is repulsive and soft, being proportional to 
1/r 3 , with r denoting the distance between the particles. The prefactor can easily be tuned 
by varying the external magnetic field strength. In real-space experiments [10,11], the two- 
stage melting process was confirmed with an intermediate hexatic phase which had a tiny 
stability range bracketed between the fluid and crystalline phase. There are also computer 
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simulations [12,13] for freezing in 1/r 3 system but finite-size effects turn out to become crucial 
here [14]. 

In this letter, we apply density functional theory (DFT) to two-dimensional freezing of soft 
1/r 3 interactions. There are two major difficulties arising in doing so: first, it is known that 
it is difficult to get a reliable density functional approximation for soft repulsive interaction 
potentials. While hard sphere freezing serves as standard test case for various approximations 
and many reliable approximations do exist (e.g. Rosenfeld's fundamental measure theory [15]), 
the freezing of soft inverse-power law-fluids turns out to be much harder [16]. For the extreme 
case of the one-component plasma, featuring a 1/r interaction, it has been shown by Likos 
and Ashcroft [17] that an extended modified weighted density approximation (EMA) which 
contains correct triplet correlation of the fluid yields reliable freezing data. In this letter, we 
overcome the first problem in a similar way and generalize the EMA to two dimensions and 
apply it to the 1/r 3 interaction. The second, more principal problem is linked to the two- 
dimensional character of the system itself. It is not clear to date how to include the hexatic 
phase into the density functional language. Here, we do not address this deep question but 
rather focus on the prediction of the freezing transition point neglecting the tiny stability 
range of the hexatic phase. A similar view has been taken for hard disk systems [18-21] and 
to the 1/r interaction in 2D [21] where density functional theory of freezing was applied to. 

We demonstrate that the EMA yields excellent freezing data as compared to the stan- 
dard modified weighted density approximation (MWDA) [22] . Even the relative- mean-square- 
displacement in the coexisting solid is in reasonable agreement with the experimental data on 
colloids in a magnetic field. 

The Helmholtz free energy density functional is typically split into the ideal gas and an 
excess part, -Ftot [p(f)] — -Fid [p{r)] + F cyL [p(r)]. Here the ideal part is local and nonlinear, 
-Fid [p(i")] = ksT J drp(r) {in [p(r)A 3 ] — l}, with A denoting the thermal de Broglie wave- 
length and ksT the thermal energy. The excess part can only be calculated approximately. 
Both the MWDA and the EMA approximate the excess free energy of the inhomogeneous 
system by setting it equal to the excess free energy of a uniform liquid evaluated at an appro- 
priately weighted density p: 

F c * [p(r)} « F C T DA / EMA W)] = Nf (p), (1) 

where N is the number of particles in the system and fo{p) is the excess free energy per 
particle of the uniform liquid at the weighted density 

/ 5 b( r )]=T7 / drdr'p(r)p(r')w(r-r';p) 

1 f &> 

+ — - / dr dr' dr"p(r)p(r')p(r")v (r — r',i — r ; p) . 
N z J 

Here the second term only appears in the EMA and not in the MWDA. The weight functions 
w(r; p) and v(r,r';p) are determined by requiring that in the uniform limit, p{r) — > p the 
approximate functional F ex [p{r)\ is exact up to second (MWDA) or third (EMA) order in 
density difference Ap(r) = p(r ) — p in the functional expansion of the excess free energy of the 
inhomogeneous system about the excess free energy of the fluid. The second weight function 
v in eq. J2J) is chosen to be zero in the MWDA since the third order correlation function does 
not enter the formalism explicitly; rather, all higher terms are approximately included as a 
consequence of the self-consistency requirement on the determination of p, appearing on both 
sides of eq. (J2J- The EMA, on the other hand, is exact up to third order and, similarly, includes 
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approximate contributions from all higher-order terms. The normalized weight functions have 
to fulfill the requirements [17] 
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where c and c are the two- and three-particle direct correlation functions of the liquid [23] 
which are an input to the theory. These conditions uniquely determine the weight functions 
and lead to simple algebraic equations for v and w that can be found in reference [17]. 

In order to find the equilibrium onc-particlc density p C q(f) we minimize the total free 
energy functional F tot [p( r )] with respect to the inhomogeneous one-particle density 

p{r). We make the Gaussian ansatz, p(r) = — J2r ex P — a \ r ~ Ri\ where {Ri} is the set 
of Bravais lattice vectors of a triangular lattice (with average density p). For fixed average 
density p we thus end up with only one minimization parameter a which describes the degree 
of localization. For a — > the density profile becomes flat and the system turns into a 
homogeneous liquid of number density p, whereas increasing a leads to enhanced particle 
localization around the lattice sites. 

With the Gaussian parametrization of the density profiles, we obtain the following self- 
consistent equation for the weighted density p as function of the localization parameter a and 
the average density p: 
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where pk = e 



(4) 
k /4a are ^ e Fourier coefficients of the Gaussian ansatz for p(r) and, likewise, 

Cq ', n = 2,3, denote the Fourier transforms of the n-particle direct correlation functions, 
evaluated at the reciprocal lattice vectors (RLV's) K and Q. Primes denote derivatives with 
respect to density and the three-particle term only appears in the EMA. Whereas Fid grows 
with a, p decreases with the latter, causing a concomitant decrease in F cx , because the latter 
is a monotonically increasing function of p. As can be induced from eq. iQJ above, the decrease 

(2) 

of p with a is pronounced when the RLV's of the lattice lie close to the maxima of Cq (fc; p), 
a feature that corresponds physically to an inherent tendency of the fluid to enhance density 
waves at these wavevectors. 

We now apply the MWDA/EMA to the inverse-power pair potential u(r) = uq/t 3 , where 
uq is a parameter with dimensions of energy x volume; for the specific realization of two- 
dimensional paramagnetic colloids of susceptibility \ exposed to a perperndicular magnetic 
field B, we have uq — (yB) 2 /2. The thermodynamics and structure depend, due to simple scal- 
ing, only on one relevant dimensionless coupling parameter T = uop 3 ^ 2 /(fcgT). Therefore, it is 
convenient to express all quantities in terms of T and consider coupling parameters rather than 
densities via this scaling relation. Correspondingly, in what follows we employ the weighted 
coupling constant T, related to p via the scaling relation T(r, a) — u p 3 ^ 2 (p, a)/(fceT). 

In order to obtain the concrete form of the functional approximations, we need the two- 
and three-particle direct correlation functions and the excess free energy per particle /o of the 
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Fig. 1 - The Fourier transform pc Q (k) of the two-particle direct correlation function at F = 9, plotted 
against ka, where a = p -1 . Shown are simulation data using the Verlet closure (solid line); liquid 
integral equation theory using the RY closure (dashed line) and liquid integral equation theory using 
the HNC closure (dotted line) . The arrows indicate the positions of the first four reciprocal lattice 
vectors of the triangular lattice. 

corresponding uniform fluid for a wide range of coupling constants I\ These quantities are 
obtained as described below: 

(i) The two-particle direct correlation function is obtained by liquid integral equation the- 
ory, where we used the hypernetted chain (HNC) [23] or the thermodynamically consistent 
Rogers- Young (RY) closure [24] . We have also produced "exact" data for the real-space direct 
two-particle correlation by computer simulation using the Verlet closure [25]. A comparison 
between the HNC, RY and simulation data for the Fourier transform Cq of the direct corre- 
lation function is shown in fig.^for the experimentally determined coupling close to freezing. 
The HNC closure underestimates the structure strongly while the RY closure is closer to the 
simulation data. We also show the positions of the first four reciprocal lattice vectors of a tri- 
angular lattice with same density. The value of c Q at these lattice vectors crucially influences 
the solid free energies, as can be seen from eqs. (QJ and Q. Therefore it can be anticipated 
that the HNC closure will not produce reliable results and we do not consider it further. 

(ii) The excess free energy per particle /o is obtained from the pair correlation via the 
compressibility route. 

(iii) Finally, the three-particle direct correlation function c n of the underlying fluid is 
obtained from an approximation by Denton and Ashcroft [26], which is based on a weighted 
density approximation to the first order direct correlation function of an inhomogeneous sys- 

(3) 

tern. This approach leads to an analytic expression of Cq in terms of the one- and two-particle 
correlation functions cj, , c Q of the liquid and their derivatives with respect to density, em- 
ploying the 'symmetrized sum' 
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Here, primes denote derivatives with respect to density. 
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Fig. 2 - The weighted coupling constant f (F, a) as a function of the localization parameter a within 
the MWDA (solid line) and within the EMA (dashed line) using the "exact" pair structure from 
simulation for the strong coupling F = 9. 

Fig. 3 - Relative free energy per particle JV -1 [i J tot(r, a) — Ftot(r, a = 0)] in units of ksT as a function 
of the localization parameter a obtained within the EMA for three different coupling constants T — 
9, 9.4, 9.8 (the three lower curves from top to bottom) compared to the relative free energy density 
obtained within the pure MWDA for a coupling constant F — 9 (uppermost line). 



Results for the approximations proposed are presented in figs. |2 and |3| In fig. |2 the 
weighted coupling constant r(F, a) is shown versus the localization parameter a for a strong 
coupling r close to freezing. Both the MWDA and the EMA are examined with the simulation 
pair structure input. Obviously, T coincides with the bare F in the fluid (a = 0). It can be 
seen that the MWDA yields a smaller reduction in T relative to F than the EMA: explicit 
inclusion of three-body effects enhances the tendency of the particles to localize. Hence, one 
expects freezing at lower couplings in the EMA. In fact, in fig. [3J the free energy difference 
between a solid of localization a and a fluid (a = 0) shown versus a reveals that the fluid 
is much more stable in the MWDA as compared to the EMA. The EMA yields a transition 
from the fluid to the solid close to F = 9.4: while for F = 9.0 the fluid is stable as indicated 
by the minimal value at a = 0, fluid-solid coexistence is achieved at F = 9.4, see the two 
equal minima in fig. [3] The solid phase, on the other hand, is clearly stable for F = 9.8. The 
localization parameter at coexistence is roughly a m [ n a 2 = 100. 

The full results of a numerical calculation using Maxwell's double tangent construction 
yield a weak first-order transition with a fluid density corresponding to a coupling constant 
of Tf and a solid density corresponding to a coupling constant of F s . There is a small density 
gap Ar = r s — Ff, describing the coexistence region. Table I summarizes the freezing/melting 
parameters for the MWDA with RY closure, for the EMA with RY closure, and for the EMA 
with the "exact" pair structure obtained from simulation. The data are compared against 
experimental results obtained from real-space microscopy measurement of magnetic colloids 
confined to an air-water interface. The experiments give freezing with an intermediate hexatic 
phase. The liquid-solid transition has also been studied using numerical simulation [12, 13] 
yielding a slightly higher inverse transition temperature between 12.0 and 12.25 but these 
investigations suffer from finite size effects. 

As becomes evident from Table I, the MWDA is not a quantitatively satisfying theory 
as it overestimates the freezing coupling by a factor of 4. Note that the overestimation of 
the freezing coupling is the reason why it is not possible to feed the "exact" pair structure 
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Table I - Freezing and melting parameters Tf and T s , the widths of the coexistence regions Ar = 
r s — T f, and the relative displacement parameters 7 at coexistence obtained within: the MWDA with 
the RY closure (first row); the EMA with the RY closure (second row); the EMA with the "exact" pair 
structure from simulation (third row) and experimental parameters for the isotropic-hexatic transition, 
the hexatic- crystal transition, and the Lindemann parameter, obtained from real-space microscopy 
measurements of magnetic colloids confined to an air-water interface (last row). 





T/ 


r s 


Ar 


7 


MWDA with RY 


41.07 


41.13 


0.06 


0.017 


EMA with RY 


23.0 


23.08 


0.09 


0.020 


EMA with simulation 


9.33 


9.49 


0.16 


0.020 


Experiment 


10.0 


10.75 


- 


0.038 



into the MWDA. At such high coupling, no fluid pair structures are available since the fluid 
spontaneously crystallizes in the simulation. The EMA, on the other hand, yields results in 
close agreement with experimental data. 

More detailed, structural information can be extracted from the localization parameter of 
the coexisting solid. For all approximations used we find localization parameters at freezing 
in the range 99 < a m i n (T f)a 2 < 115. Strictly speaking, the localization parameter has no 
counterpart in "real" 2D systems since the particles are not localized due to long range fluctu- 
ations. However, if one relates the particle displacements to that of their nearest neighbor, one 

can define a finite quantity as 7 = p < (Ui — Uj_|_i) V where Ui and itj+i are the displacement 

vectors of neighboring lattice sites. Disregarding nearest-neighbor correlations (ut ■ Uj+i), 
7 can be estimated. Since the nearest- neighbor correlations (v,i -Wj+i) are expected to be 
positive: 

7<2 /0 ( M 2)«2/(a min a 2 ). (7) 

By this relation, the localization parameter of the coexisting solid gives a prediction for 7 which 
is included in Table I. From experiments, 7 is known to be close to = 0.038 [10]. This was 
shown to be in accordance with harmonic lattice theory [27]. The EMA yields 7 < 0.020, i.e. 
the EMA roughly overestimates the localization of the particles by a factor of 2. 7 is smaller 
than the experimental value, contrarily to what was expected from the inequality 0. This 
shows that there is still a need to improve the theories in order to correctly predict localization 
properties. A similar overestimation of the localization is also common in weighted density 
approximations in three spatial dimensions [22]. 

In conclusion, we have demonstrated that the EMA is able to quantitatively predict the 
freezing transition of a two-dimensional colloidal system with soft and long-ranged 1/r 3 - 
intcractions in good agreement with experimental and simulation data. In analogy to three- 
dimensional systems, the appearance of long-range interactions requires the explicit inclusion 
of three-particle correlation functions of the liquid in the construction of the weighted density. 
Furthermore, the predicted transition temperatures are very sensitive towards slight changes 
of the two-particle correlation functions of the underlying fluid. A highly accurate input of 
the same is therefore crucial. 

Relying on the good quality of the EMA functional, our results can serve as a platform to 
treat more challenging problems than bulk transitions [28] . One obvious extension is towards 
external potentials acting on the particles, such as system walls or gravity. The MWDA 
can in principle be applied to a fluid near a single wall, but not to a free interface between 
coexisting phases. Another interesting example is a spatially inhomogeneous magnetic field, 
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which renders the interactions space dependent [29]. Finally, one may employ the scheme of 
dynamical density functional theory [30, 31] in order to use the EMA functional to study the 
effect of spatially homogeneous magnetic fields that oscillate in time. Work along these lines 
is currently under way. 
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